\documentclass[11pt]{article} 
\usepackage{amsmath,amsthm,amstext,amssymb,amsfonts,amscd,graphicx,bbm,algorithmic} 

\voffset -0.5in
\addtolength{\textheight}{1in} 

% new theorems
\newtheorem{lem}{Lemma} 
\newtheorem{lemma}{Lemma}
\newtheorem{thm}{Theorem} 
\newtheorem{ass}{Assumption} 
\newtheorem{cor}{Corollary} 
\newtheorem{clm}{Claim}
\newtheorem{prop}{Proposition} 
\newtheorem{con}{Conjecture} 
\newtheorem{wff}{Wff} 
\newtheorem{defn}{Definition} 
\newtheorem{axiom}{Axiom} 
\newcounter{rulenum} 
\newcounter{sent} 
\newcounter{tempcnt} 

% macros
\newcommand*{\add}{\mbox{\bf add}} 
\newcommand*{\open}{\mbox{\bf open}} 
\newcommand*{\close}{\mbox{\bf close}} 
\newcommand*{\ex}{\mbox{\rm E}} 
\newcommand*{\pr}{\mbox{\rm Pr}} 
\newcommand*{\range}{\mbox{\rm range}} 
\newcommand*{\rank}{\mbox{\rm rank}} 
\newcommand*{\sgn}{\mbox{\rm sign}} 
\newcommand*{\var}{\mbox{\rm Var}} 
\newcommand*{\diag}{\mbox{\rm diag}} 
\newcommand*{\epi}{\epsilon} 
\newcommand*{\QED}{\ \hfill\rule[-2pt]{6pt}{12pt} \medskip}
\newcommand*{\supp}{\mbox{\rm supp}} 

\newcommand*{\grad}{\nabla}
\newcommand*{\half}{\frac{1}{2}}
\newcommand*{\inv}{^{-1}}
\newcommand*{\0}{\mathbf{0}}
\newcommand*{\1}{\mathbf{1}}
\newcommand*{\E}{\ensuremath{\operatorname{E}}}
\newcommand*{\maximize}{\text{maximize}}
\newcommand*{\minimize}{\text{minimize}}
\newcommand*{\st}{\text{subject to}}
\newcommand*{\R}{\mathbbm{R}}
\newcommand*{\matlab}{{\sc Matlab}}

\newcommand{\abs}[1]{\left\vert #1 \right\vert}
\newcommand{\bigo}[1]{\mathcal{O} \left( #1 \right)}
\newcommand{\cov}[2]{\ensuremath{\operatorname{Cov}\left( #1, #2\right)}}
\newcommand{\Ex}[2][]{\ensuremath{\E_{#1} \left[ #2 \right]}}
\newcommand{\norm}[1]{\left\lVert\,#1\,\right\rVert}
\newcommand{\bmat}[1]{\begin{bmatrix}#1\end{bmatrix}}
\newcommand{\pmat}[1]{\begin{pmatrix}#1\end{pmatrix}}
\newcommand{\smallmat}[1]{\left (\begin{smallmatrix}#1\end{smallmatrix} \right)}
\newcommand{\vb}[1]{\mathbf{#1}}

\renewcommand{\P}{\ensuremath{\operatorname{P}}}
\renewcommand{\Pr}[2][]{\ensuremath{\P_{#1} \left \{ #2 \right \}}}


\title{Existence of Positive Steady States for Mass Conserving and Mass-Action Chemical Reaction Networks with a Single Terminal-Linkage Class%
   \thanks{\today.}}
\author{Santiago Akle%
   \thanks{Institute for Computational and Mathematical Engineering,
           Stanford University, Stanford, CA 94305.
           Research supported in part by the U.S. Department of Energy (Office of
Advanced Scientific Computing Research and Office of Biological and
Environmental Research) as part of the Scientific Discovery Through
Advanced Computing program, grant DE-SC0002009.
           Email: {\tt akle@stanford.edu}, {\tt onkar@stanford.edu}, 
                  {\tt ntaheri@stanford.edu}.}
   \and Onkar Dalal\footnotemark[2]
   \and Ronan M.T. Fleming\thanks{Center for Systems Biology, University of Iceland, Sturlugata 8,
Reykjavik 101, Iceland. Email: {\tt ronan.mt.fleming@gmail.com}}
   \and Michael Saunders\footnotemark[4] 
   \and Nicole Taheri\footnotemark[2]
   \and Yinyu Ye%
   \thanks{Department of Management Science and Engineering,
           Stanford University, Stanford, CA 94305. Research supported in part by NSF grant GOALI 0800151
           and DOE grant DE-SC0002009.
           Email: {\tt saunders@stanford.edu}, {\tt yinyu-ye@stanford.edu}.}}
\date{}
\begin{document}
\maketitle

\begin{abstract}
We establish that mass conserving single terminal-linkage networks of chemical reactions admit positive steady states regardless of the choice of reaction rate constants and network deficiency. This result holds for closed systems without material exchange across the boundary as well as for open systems with material exchange at rate that satisfy a simple sufficient and necessary condition. 

Our proof uses a fixed-point of a novel convex-optimization formulation for analyzing the steady state behavior of networks of chemical reactions satisfying the law of mass-action kinetics. We suggest a fixed-point algorithm to compute these steady states and report some numerical experiments. 
\end{abstract}

\section{Introduction} 

Deterministically modeled chemical reaction networks are of interest in the field of systems biology and elsewhere. Models based on the assumption of mass-action kinetics depend on the kinetic parameters of the system. However, measuring these parameters empirically is difficult and error-prone. Hence we focus on establishing properties that are independent of the particular choice of parameters.

In this work we address the issue of existence of positive steady states,
i.e.,\ concentrations of species that will stay constant under the system's
dynamics. The particular case we tackle is where the directed graph of chemical
reactions forms a \emph{strongly connected component} and the reactions
conserve mass. We prove the existence of strictly positive equilibria for
closed systems having no exchange across the boundary and extend our methods to
open systems exchanging species at certain rates. How best to compute such a
steady state remains uncertain; however, for closed systems we suggest a
fixed-point algorithm and provide results of our numerical experiments for
several cases. 

\subsection{Background}
A mathematical theory for this problem, also known as \textit{Chemical Reaction
Network Theory} (CRNT), has its roots in the seminal work by Fritz Horn, Roy
Jackson, and Martin Feinberg in \cite{GMAK, uniqueEPandLyapunov, necc-suff-CB,
deficiency0, deficiency1}. We build on the notation used in these works and
summarized by Gunawardena et al. \cite{gunawardena}. The system consists of a
collection of species reacting collectively in some combination to give another
combination of species in a network of chemical reactions.  Let $\mathcal{S}$
be a set of $m$ species and $\mathcal{C}$ be a set of $n$ complexes. The
relation between species and complexes can be written as a non-negative matrix
$Y\in\R^{m\times n}$, where column $y_j$ represents complex $j$ and $Y_{i j}$
is the multiplicity of species $i$ in complex $j$.  For example, the
multiplicity of the species NaCl in the complex (H$_2$O + 2NaCl) is 2.  

A reaction network is represented by the underlying weighted directed graph
$G(V,E)$ with each node in $V$ denoting a complex, each directed edge
$i\rightarrow j$ denoting a reaction using $i$ to generate $j$, and the
positive weight $k_{i\rightarrow j}$ being the reaction rate.  The matrix $A
\in \R^{n \times n}$ is the weighted adjacency matrix of the graph, where
$A_{ij}=k_{i\rightarrow j}$.  Define $D := \diag(A\1)$, where $\1$ is the
vector of all ones, and $A_k := A^T-D$. The vector for concentration of species
in a reaction network is given by $c\in\R^m$, and the vector of rates of
exchange of species across the network boundary is given by $b$. Our aim is to
investigate necessary and sufficient conditions on $b$ to prove the existence
of $c$ and be able to compute it given a specific $b$. A nonlinear function
$\psi(c):\R^m_+\rightarrow \R^n$ capturing the mass-action kinetics is given by
\[\psi_j(c) = \prod_i\,c_i^{Y_{ij}},\] where the system of ordinary
differential equations governing the reaction network can be written as
\[\dot{c} = YA_k\psi(c) - b.\] \noindent Hence, a steady state concentration
for a chemical reaction network is any non-negative vector $c^*\in\R^m$ such
that $YA_k\psi(c^*)=b$. Equivalently, a pair $(c^*,v^*)$ with $c^*\in\R^m$ and
$v^*\in\R^n$ will be a steady state if it satisfies the conditions
\begin{align} 
  YA_kv^* &=b \label{fb}\tag{FB} \\ \psi(c^*) &= v^*
  \label{mak}\tag{MA} 
\end{align}
or \emph{flux-balance} and \emph{mass-action}, respectively. We can find a steady state $c^*$ by finding a vector $v^*$ that satisfies $\eqref{fb}$ and also satisfies \eqref{mak} for some vector $c^*\in\R^m$.

Observe that if $v^*$ is a positive steady state, then 
\begin{align}
  Y^T\log(c^*)= \log(\psi(c^*))= \log(v^*).
  \label{mak-alt}
\end{align} 
We refer to this alternative condition as the logarithmic form of the mass-action condition \eqref{mak}.

In this notation, systems with the property of mass conservation can be characterized by the following definition. 
\begin{defn}
A chemical reaction network is $\emph{mass conserving}$ if and only if there exists a strictly positive vector $e\in\R^m$ such that 
\begin{align}
 e^TYA_k=0,
  \label{consis}
\end{align}
where $e$ denotes the molecular weights of the species (or atomic weights if the species are elements).  
\end{defn}

The connectedness of the networks is captured in the following definition of a terminal-linkage class.  
\begin{defn} A \emph{terminal-linkage class} is defined as a set of complexes $\mathcal{L}$ such that for any pair of complexes $(i,j) \in \mathcal{L}$ there exists a directed path in the graph $G$ that leads from $i$ to $j$.  
\end{defn} 

We further restrict our analysis to a class of weakly reversible networks.
\begin{defn} A chemical reaction network is \emph{weakly reversible} if it is formed exclusively by one or more terminal-linkage classes.
\end{defn} 

A reaction network that consists of exactly one terminal-linkage class is called a \emph{single terminal-linkage network}.

As suggested by simple examples like a single non-reversible reaction, reversibility at least in a weak sense is a prerequisite for steady states with strictly positive concentrations for all species.

Another useful definition for a network is a stoichiometric subspace.
\begin{defn} A \emph{stoichiometric subspace} $S$ is the subspace defined by the span of vectors $y_{j}-y_{i}$, where $y_{j}, y_{i}$ are the columns of $Y$ denoting the complexes $i$ and $j$ for each reaction pair $i \rightarrow j$ in the network. 
\end{defn} 
In one of the early works, Horn and Jackson \cite{GMAK} analyzed mass-action
kinetics for closed systems (with $b=0$) and defined a class of equilibrium
points, called \emph{complex-balanced equilibria}, and systems admitting such
an equilibrium, \emph{complex-balanced systems}. These closed systems are shown
to satisfy the \emph{quasi-thermostatic} and \emph{quasi-thermodynamic}
conditions regardless of the kinetic rate constants. Following this, Horn
\cite{necc-suff-CB} also proved necessary and sufficient conditions for
existence of a \emph{complex-balanced} equilibrium. In
\cite{uniqueEPandLyapunov}, Feinberg and Horn used the existence of a Lyapunov
function to show the uniqueness of the strictly positive steady state in each
stoichiometric compatibility class, which is equivalent to specifying all the
conserved quantities of a system. Later Feinberg \cite{deficiency0,
deficiency1} proved two theorems, now famously known as Deficiency 0-1
theorems, that provide the analysis of strictly positive steady states for a
class of networks with deficiency $0$ or networks with deficiency $1$ but with
each terminal-linkage class having deficiency less than $1$.  The
\emph{deficiency} of a network is defined as $\delta = n - t - s$, where $t$ is
the number of terminal-linkage classes and $s$ is the dimension of the
stoichiometric subspace, also known as stoichiometric compatibility class. For
this restricted class of closed systems, the existence of a positive steady
state is given by Perron-Frobenius theory for a strictly positive eigenvector.
Other work in this area is from the perspective of dynamical systems and aimed
toward proving two open conjectures: \emph{Global Attractor Conjecture} and
\emph{Persistence Conjecture} \cite{Anderson-GAC}. Another approach using
parametrized convex optimization to compute a non-equilibrium steady state is
given in \cite{fleming-opt}.

To the best of our knowledge, the vast majority of CRNT research studies closed
\textit{complex-balanced systems}, which, by definition, are assumed to admit a
\textit{complex-balanced equilibrium}. In the notation above, complex-balanced
equilibrium holds when the vector of concentrations satisfies $A_k\psi(c)=0$,
i.e., the vector $\psi(c)$ belongs to the null space of $A_k$. It can be shown
that a network with linearly independent complexes will have deficiency $\delta
= 0$. However, if some of the complexes are linearly dependent (as shown in the
example in Section 3), there are systems that are not complex-balanced but
admit equilibrium points with vectors $\psi(c)$ such that $A_k\psi(c) \neq 0$
but $A_k\psi(c)$ belongs to the null space of $Y$. Though the condition of
\emph{complex-balance} is sufficient for \emph{thermodynamic} consistency,
\cite{GMAK} shows that it is not necessary. Also, for open systems, those with
material exchange across their boundary,\emph{complex-balance} is not defined.
In order to handle open systems, these works hint at extending the system using
a \emph{pseudo $0$-complex} and adding \emph{pseudo reactions}. However, it is
unclear how to choose the \emph{pseudo kinetic rates} such that the positive
eigenvector solution of the extended system will achieve the given external
exchange rates $b$. From the point of view of systems-biology and bio-chemical
engineering, analyzing the behavior of a cell under different exchange
conditions $b$ is very important to control and engineer the cell for e.g.
studying desired effects in pharmacology or producing specific metabolites in
bioreactors.

\vspace{0.2in}
In this paper, we extend the previous work on two accounts: 1) extending the
proof for existence of positive equilibria for closed systems to include some
systems that do not satisfy the necessary conditions of the Deficiency 0-1
theorems (are not necessarily \emph{complex-balanced}), and 2) proving a
necessary and sufficient condition on the external exchange rate $b$ for some open
systems to admit a strictly positive steady-state. We use a fixed-point of a
mapping based on a convex optimization problem. The objective function of the
optimization problem is very close to the Helmholtz function defined in
\cite{GMAK}. The fixed point of this mapping gives the required steady state
solution. We prove the existence of a strictly positive steady state for any
weakly reversible chemical reaction network with a \emph{single
terminal-linkage class}. We strongly believe that this can be extended to
systems with \emph{multiple terminal-linkage classes}, as supported by our
computational results for randomly generated networks. Section 3.2 gives a
detailed analysis of a toy network to emphasize our extension.


\section{A Fixed Point Model} 

Our main result establishes that for any set of reaction rates $k>0$ and 
any $b$ in the range of $YA_k$, a network formed by a single terminal-linkage
class will admit a strictly positive solution pair $(c,v)$ that satisfies the
laws \eqref{fb} and \eqref{mak}.  We show this by defining a related mapping
with positive fixed points, and establishing an equivalence between these
positive fixed points and positive solutions to the equations.

We define a mapping in terms of a linearly constrained optimization problem.
The problem is constructed so that the logarithmic form of the mass-action
equation \eqref{mak-alt} is an optimality condition, and hence any solution to
this optimization problem will also satisfy \eqref{mak-alt}.

To define the mapping, let $b=YA_k\eta$ for some $\eta$, and observe we can
write $b = YA^T(\eta +s) - YD(\eta + D^{-1}A^Ts)$, with $s$ chosen arbitrarily.
In particular we can choose $s$ positive and large enough so that $\eta^+ :=
\eta+s$ and $\eta^- := \eta + D^{-1}A^Ts$, are both positive. Also from 
Definition \eqref{consis} $e^Tb = 0$ thus
\begin{equation}
  e^TYD\eta^- = e^TYA^T\eta^+.
  \label{massbalanceb}
\end{equation}
Let $\mu:=(r,r_0) \in \R^{m+1}$ be a vector parameter.  Observe that if the
parametric convex optimization problem
\begin{align}
  \underset{(v,v_0)\in\R^{n+1}}{\minimize} &\quad v^TD(\log(v)-\1) + v_0(\log v_0 -1)   \notag
\\                     \st &\quad YD v + YA^T\eta^+ v_0 = YA^Tr + YD\eta^-r_0 &:\ y \label{convex-fix}
\\                         &\quad (v,v_0) \ge 0                     \notag
\end{align}
has a strictly positive solution $(v^\star(\mu),v_0^\star(\mu))$, then the
optimality conditions
\begin{align}
   YDv^\star(\mu)  + YA^T\eta^+ v^\star_0  &= YA^Tr + YD\eta^-r_0  \notag
\\ DY^T y^\star(\mu) &= D\log(v^\star(\mu))  \label{optcon}
\\ (YA^T\eta^+)^Ty^\star(\mu)   &= \log(v^\star_0(\mu)) 										  \notag
\\                   (v^\star(\mu),v_0^\star(\mu)) &\ge 0                    \notag
\end{align} 
are well defined. Since $D$ is nonsingular, the second optimality condition is
equivalent to \eqref{mak-alt}, and hence implies mass-action. We show that
for some parameter $\mu$ the equality 
\begin{equation}
 YA_kv^\star=v^\star_0b
  \label{scaled-fb}
\end{equation} is
satisfied. This implies that for $b=0$ the solution is attained. For the case
$b\neq 0$, we show that we construct a second solution so that $\eqref{fb}$ is
attained. 

Note that the nonlinear program \eqref{convex-fix} is strictly convex, so for any
feasible $(r,r_0)$ there is a unique minimizer. That is, the mapping 
\begin{equation}
  (r,r_0) \rightarrow (v^\star(\mu),v^\star_0(\mu))
  \label{mapping}
\end{equation} is well defined. 

If $\mu$ is a fixed point of \eqref{mapping}, then the linear
equality constraint in \eqref{convex-fix} implies
\[
   YDv^\star(\mu)+YA^T\eta^+v^\star_0(\mu) = YA^Tv^\star(\mu) + YD\eta^-v_0^\star(\mu)
\]
or, equivalently,
\[
    Y A_k v^\star(\mu) = Y(A^T-D)v^\star(\mu) = v^\star_0(YA^T\eta^+ - YD\eta^-) = v^\star_0(\mu)b
\]
Therefore, if such a fixed point exists, the solution $v^\star(\mu)$ at this
fixed point will satisfy $\eqref{scaled-fb}$.  For simplicity, we henceforth refer to
the optimal solution variables $(v^\star(\mu),v^\star_0(\mu)), y^\star(\mu)$ as
$(v^\star,v^\star_0), y^\star$, but acknowledge their dependence on $\mu$.

\begin{thm} \label{fp-exist-map} For any mass conserving, mass-action chemical
  reaction network and any choice of rate constants $k>0$, there exist
  nontrivial fixed points for the mapping given by \eqref{mapping}.
\end{thm} 

\begin{proof} 
	Brouwer's fixed point theorem states that any continuous mapping from a
	convex and compact subset of a Euclidean space $\Omega$ to itself must have
	fixed points. 

	Let $(v^\star,v^\star_0)$ be defined as in \eqref{mapping} and let
	$\gamma>0$ be a positive and fixed scalar. Define the set \[\Omega =
	\left\{ (v,v_0)\in \Re^{n+1} \; : \; (v,v_0)\geq 0, \quad e^TYDv +
	e^TYA^T\eta^+v_0 = \gamma   \right\}.\]  where $e$ is defined in
	\eqref{consis}.

	The set $\Omega$ is formed by an intersection of closed convex sets and is
	bounded, therefore it is convex and compact. The mapping $\mu\rightarrow
	(v^\star,v^\star_0)$ is continuous. Since problem \eqref{convex-fix} is
	feasible for any $\mu\in\Omega$, the mapping $\Omega \ni \mu \rightarrow
	(v^\star,v^\star_0)$ is well defined.

	To show the image of $\Omega$ under the mapping is in $\Omega$, first
	observe that by the bounds in \eqref{convex-fix}, $(v^\star,v^\star_0) \ge 0$.
	Using the equality constraints, Definition \eqref{consis} and equation \eqref{massbalanceb}
	\begin{align}
	e^T YD v^\star + e^T YA^T \eta^+ v^\star_0 &= e^T YA^T r + e^T YD \eta^- r_0 \notag 
	\\ & = e^T YD r + e^T YA^T \eta^+r_0=\gamma, \notag \end{align} thus
	$(v^\star,v_0^\star)\in \Omega.$

	To conclude, the mapping $\Omega \ni \mu \rightarrow (v^\star,v^\star_0) \in
	\Omega$ exists and must have a fixed point.  Moreover, since $\Omega$ does
	not contain the zero vector, the fixed point(s) are nontrivial.
 
  \end{proof}

We have established the existence of a nontrivial fixed point $\mu$ of the
mapping $\Omega \ni \mu \rightarrow (v^\star,v^\star_0)\in \Omega$. We have also
argued that when the associated minimizer $(v^\star,v^\star_0)$ is strictly
positive, it is a solution to \eqref{mak} and to $YA_kv^\star=v^\star_0b$.
However, in the case when some entries of $v^\star$ are zero, the objective
function of \eqref{convex-fix} is non-differentiable and we cannot use the
optimality conditions to show that \eqref{mak} holds. 

The value of $YDr+YA^T\eta^+r_0$ is the rate of consumption of each chemical
species and $YA^Tr+YD\eta^-r_0$ is the rate of production of each chemical
species. At the fixed point, the equality $YDr + YA^T\eta^+r_0=
YA^Tr+YD\eta^-r_0$ defines a steady state. Observe that the definition of
$\Omega$ involves the parameter $\gamma=e^T(YDr + YA^T\eta^+r_0)$.  Since the
vector $e$ can be interpreted as an assignment of relative mass to the species,
$\gamma$ can be interpreted as the total amount of mass that reacts per unit
time at the steady state.  Therefore looking for fixed points in $\Omega$
corresponds to looking for steady states where the ammount of mass that reacts
in the system is prescribed.
	

\subsection{Positive fixed points in single terminal-linkage networks}  \label{section::single-linkage}

We now consider the case when the network is formed by a single
terminal-linkage class and show that if $\hat \mu$ is a fixed point of the
mapping \eqref{mapping}, the minimizer $(v^\star(\hat
\mu),v^\star_0(\hat\mu))$, and therefore $\hat\mu$, is strictly positive.

Lemma $\ref{maximum-support}$ shows that if problem $\eqref{convex-fix}$ has a
feasible point with support $J$, the minimizer $(v^\star,v^\star_0)$ will have
support at least $J$. Lemma \ref{positive-feasible} expands on this result and
uses the single terminal-linkage hypothesis to show that at a fixed point,
there is a positive feasible point.  These two Lemmas imply that at a fixed
point $\hat{\mu}$, the minimizer will be strictly positive.
Since the minimizer is the fixed point, the fixed point will be positive.
Finally Theorem \ref{thm:scaling} shows that if at the solution $\hat{v}_0 \neq
1$ we can construct another solution for which $\hat{v}_0=1$. This establishes
that there is a nontrivial steady state for the network. 

To complete the argument we must establish Lemmas \ref{maximum-support},
\ref{positive-feasible} and Theorem \ref{thm:scaling}.

\begin{lemma} \label{maximum-support}
  The minimizer $(v^\star,v^\star_0)$ of Problem \eqref{convex-fix} 
  has support at least as large as any feasible point.
\end{lemma}
\begin{proof}
  Let $\tilde{v}\in \Re^{n+1}$ be any of the feasible points with the largest
  support and let $z$ be any feasible direction at $\tilde{v}$. By
  construction, for all $\alpha$ in some interval $[\ell,u]$ the points
  $v_\alpha:= \tilde{v} + \alpha z$ are feasible. Notice that WLOG we can
  assume that $\ell<0$ and $0<u$ since if not $\ell = 0$ and $u>0$, in which
  case we can choose $v_{u/2}$. Since $v_{u/2}$ is a convex combination
  of $\tilde{v}$ and $\tilde{v}+uz$, it has support as large as $\tilde{v}$.

  The interval can be chosen so that when $\alpha=l$ and when $\alpha=u$, one new
  bound constraint becomes active.  This implies that $\supp(v_\ell)$ and
  $\supp(v_u)$ are strictly contained in $\supp(\tilde{v})$.

	Define the univariate function \begin{equation}   \label{univariate}
	  g(\alpha) := \phi(\tilde{v} + \alpha z), \end{equation} where $\phi$ is
	  the objective function of \eqref{convex-fix}.  We will establish that as
	  $\alpha\rightarrow l$ the derivative $g'(\alpha)\rightarrow -\infty$, and
	  as $\alpha\rightarrow u$ the derivative $g'(\alpha)\rightarrow \infty$.
	  Thus, by the mean value theorem, there must exist a zero in the interior
	  of the interval $[l,u]$.  Since this function is strictly convex, if a
	  stationary point exists in the interior of the interval, the function
	  value at the stationary point must be smaller than at the boundary.

  Observe that if we let $d_i$, for $i\in[1,\dots n]$ be the diagonal entries of $D$ and 
  $d_{n+1}=1$, we can write 
  \begin{align*}
	 g(\alpha) &=\sum_i (\tilde{v} + \alpha z_i) d_i\log(\tilde{v}_i + \alpha z_i).
   \end{align*} 
   An important observation is that if some entry $\tilde{v}_j=0$ then $z_j=0$,
   otherwise $v_\alpha$ would have a larger support for some $\alpha$. This
   implies that for all such entries $(v_\alpha)_j=0$.  If we let $J$ be the
   set of non zero entries of $\tilde{v}$ then we can write   
   \begin{align*} g'(\alpha) &=\sum_{i\in J} z_i d_i(\log(\tilde{v}_i + \alpha
	 z_i)) \\            &= \sum_{i\in L^c\cap J}
	 z_id_i(\log{(\tilde{v_i}+\alpha z_i)}) + \sum_{i\in L}
	 z_id_i(\log{(\tilde{v_i}+\alpha z_i)}), \end{align*} where $L$ is the
	 subset of $J$ formed by the entries that tend to zero as $\alpha
	 \rightarrow l$. As $\alpha \rightarrow l$, the first summation will
	 approach a finite value.  Since $z_i>0$ for all $i\in L$, the entries in
	 the logarithm of the second sum tend to zero and the term will diverge to
	 $-\infty$.

	Let $U$ be the subset of $J$ formed by the entries which tend to zero as 
	$\alpha\rightarrow u$. Observe that for these entries $z_i<0$. We can
	write
  \[
    g'(\alpha) = \sum_{i\in U^c\cap J} z_i d_i(\log{(\tilde{v_i}+\alpha z_i)})
               + \sum_{i\in U}   z_i d_i(\log{(\tilde{v_i}+\alpha z_i)}),
  \]
	in which the first sum will tend to a finite value and the second will
	diverge to $\infty$. Therefore, there must exist a stationary point strictly
	in the interior of the interval $[l,u]$.

	Now, assume that for some $\mu$ there is a feasible point $(\tilde{v},\tilde{v_0})$ with
	larger support than the minimizer $(v^\star,v^\star_0)$ of problem \eqref{convex-fix}.
	By construction, the direction $z = (v^\star,v_0^\star)-(\tilde{v},\tilde{v_0})$ is feasible, and by
	the previous argument there is a point in the line spanned by $z$
	containing $(\tilde{v},\tilde{v_0})$ that has a lower function value than $(v^\star,v^\star_0)$,
	contradicting the minimality of $(v^\star,v^\star_0)$.
  \end{proof}



\begin{lemma}\label{positive-feasible}
  If the network is formed by a single terminal-linkage class, when Problem \eqref{convex-fix} is parametrized by the fixed point $\hat\mu$,  there exists a strictly
  positive feasible point $(\hat v,\hat v_0)$.
\end{lemma}

\begin{proof}
  Let Problem \eqref{convex-fix} be parametrized with the fixed
  point $\hat\mu$, and let $(\hat v, \hat v_0)$ be both the minimizer and the
  fixed point.  We prove by contradiction that no entry of the minimizer
  $(\hat{v},\hat{v}_0)$ can be zero. Observe that by the definition of $\Omega$
  the origin is not contained in the set, therefore the fixed point cannot be identically zero. 

	Assume that $\hat{v}_0>0$, and observe that $\rho := (D^{-1}A^T\hat v + \eta^-\hat v_0, 0)$ is 
	a feasible point. Since $\eta^-$ was chosen to be strictly positive and 
	$D^{-1}A^T$ has no zero columns, the support of $\rho$ are the first $n$ entries 
	of the vector. A convex combination of $(\hat v,\hat v_0)$ and $\rho$ will have 
	full support, and will be feasible. 

	Assume that some entry of $\hat v$ is non zero and that $\hat v_0=0$, and
	observe that $(D^{-1}A^T\hat{v},0)$ is feasible. A convex combination
	of $(D^{-1}A^T\hat{v},0)$ and $(\hat{v},\hat{v_0})$ is feasible and its
	support contains the union of the supports of the two vectors. The support
	of $D^{-1}A^T\hat v$ is the support of $A^T\hat v$. Therefore by Lemma
	\ref{maximum-support},
	\begin{equation} \supp(\hat{v},\hat{v_0}) \supset \supp(\hat{v},\hat{v}_0) \cup \supp(A^T\hat{v},\hat v_0).
	  \label{monotonesets} 
	\end{equation} This relation can be used inductively
	  to show that the support of $\hat{v}$ is at least as large as the union
	  of the supports of $((A^T)^p\hat{v})$ for all positive powers of $p$.

	The single terminal-linkage class hypothesis implies that for any pair
	$(i,j)\in [1,\dots,n]\times [1,\dots,n]$, there exists a power $p$ such that
	$(A^T)^p_{ij}>0$.  Thus, if $j$ is in the support of $\hat v$, then $i$ is in
	the support of $(A^T)^p\hat v$ and hence $i$ is in the support of $\hat v$.
	Therefore, if $\hat v \neq 0$, then $\hat v>0$.

	Finally if $\hat v > 0$ there is a scalar $0<\alpha$ small enough such that 
	\[0< \hat v - D^{-1}A^T\eta^+ \alpha.\]
	Therefore the equality  
	\[
		YD(\hat v - D^{-1}A^T\eta^+\alpha)+YA^T\eta^+\alpha = YD\hat v = YA^T\hat v
	\]
    implies that the strictly positive point $(\hat v - D^{-1}A^T\eta^+\alpha,\alpha)$ is feasible.
\end{proof}

\begin{thm}\label{thm:scaling}
  If the network is mass conserving and formed by a single terminal-linkage
  class, then for any $b$ in the range of $YA_k$ there exists a positive
  concentration $c>0$ such that $YA_k\psi(c) = b$. 
\end{thm}
\begin{proof}
  Assume that there exist a positive $c>0,\alpha>0$ such that $YA_k\psi(c) =
  \alpha b$.  If we can construct $\tilde c >0$ such that $\psi(\tilde c) =
  \frac{1}{\alpha}\psi(c)$, then \[YA_k\psi(\tilde c) =
  \frac{1}{\alpha}YA_k\psi(c) = b.\]

  Let $\1$ denote the vector of all ones and observe that $\log\left(
  \frac{1}{\alpha}\psi(c) \right) = \log(\psi(c))-\1\log(\alpha) = Y^T\log(c) -
  \1\log(\alpha)$. Assume that $\1$ is in the range of $Y^T$ and that can be
  expressed as $Y^T\delta = \1$.  Then $\log(\alpha)\1 =
  Y^T(\log(\alpha)\delta)$ thus if we let $\log(\tilde c) = \log(c) -
  \log(\alpha)\delta$ (here $\log(c)$ and $\log(\tilde c) $) are entry-wise
  logarithms of the corresponding vectors and $\log(\alpha)$ the scalar
  logarithm), then \[Y^T\log(\tilde c) = Y^T\log(c) - \log(\alpha)\1,\] which
  implies that \[\psi(\tilde c) = \frac{1}{\alpha}\psi(c),\] and thus the
  inhomogeneous system can be solved.

  Finally we have to argue that the vector $\1$ is in the range of $Y^T$ for
  single terminal-linkage-class networks. Observe that the condition of mass
  conservation implies that $e^TY A_k = 0$. This implies that $Y^Te\in
  {\mathcal N}(A_k^T)$.  Since the matrix $A_k^T$ is the graph laplacian of a
  strongly connected graph it contains only the constant vector in its null
  space, thus $Y^Te = \beta\1$ for some $\beta$.  \end{proof}
 
\input{section3}
\frenchspacing
\nocite{*} \bibliographystyle{plain} \bibliography{fixpoint}{}

\end{document}
